'''
Created on Dec 1, 2010

@author: mkiyer

This module contains classes and parsers for pipeline configuration parameters
stored in configuration files.  Currently, configuration files are in XML format
''' 
import os
import re
import collections
import logging
import xml.etree.cElementTree as etree
import ConfigParser

from base import parse_bool

class PipelineConfig(object):
    '''
    this is the base class for pipeline configuration.  different
    pipelines should inherit this class and add pipeline-specific
    parsing
    
    the default pipeline provides the capability to:
    (1) copy files,
    (2) uncompress files
    (3) run the FastQC profiling tool
    '''    
    def from_xml(self, xmlfile):
        tree = etree.parse(xmlfile)        
        root = tree.getroot()
        # directories
        self.output_dir = root.findtext('output_dir')
        self.tmp_dir = root.findtext('tmp_dir') 
        # fastqc tools
        self.fastqc_bin = root.findtext('fastqc_bin')

class ExomePipelineConfig(PipelineConfig):
    '''
    exome-specific pipeline configuration
    '''
    
    def from_xml(self, xmlfile):
        # parent class init
        super(ExomePipelineConfig, self).from_xml(xmlfile)
        tree = etree.parse(xmlfile)    
        root = tree.getroot()
        
        # subdirectories for different components of the pipeline
        self.analysis_dir = os.path.join(self.output_dir, "analysis")
        self.sample_dir = os.path.join(self.output_dir, "sample")
        self.lane_dir = os.path.join(self.output_dir, "lane")
        
        # email addresses
        self.email_addresses = []
        for elem in root.findall("email"):
            self.email_addresses.append(elem.text)
            
        # Run parameters
        
        self.run_param = {}
        for elem in root.findall("run"):
            thisrun = elem.get("name")
            do_recal = elem.findtext('do_recalibration')
            quick_dirty = elem.findtext('quick_dirty')
            paired_samples =  elem.findtext('paired_samples')
            multi_samples = elem.findtext('multi_samples')
            analysis_type = elem.findtext('analysis_type')
            self.run_param[thisrun] = type('RunParam',(object,),
                                           dict(name=thisrun,
                                                do_recal=eval(do_recal),
                                                quick_dirty=eval(quick_dirty),
                                                paired_samples=eval(paired_samples),
                                                multi_samples=eval(multi_samples),
                                                analysis_type=analysis_type,
                                                )
                                           )
            
        # tool paths
        #self.samtools_path = root.findtext('samtools_path')
        self.picard_path = root.findtext('picard_path')
        self.rscript_bin = root.findtext('rscript_bin')
        self.vcftools_path = root.findtext('vcftools_path')
        self.path_to_fastx = root.findtext('fastx_path')
        
        # bwa params
        elem = root.find('bwa')
        self.bwa_path = elem.get('path')
        self.bwa_use_mem = float(elem.findtext('use_mem'))
        self.bwa_num_cores = int(elem.findtext('num_cores'))
        self.bwa_num_mismatches = int(elem.findtext('num_mismatches'))
        
        # sam tools params
        elem = root.find('sam')
        self.samtools_path = elem.get('path')
        sam_min_mapping_qual=elem.findtext('min_mapping_quality')
        sam_min_base_qual=elem.findtext('min_base_quality')
        self.sam_mpileup_params=[sam_min_mapping_qual,sam_min_base_qual]
        
        self.sam_target_exons = elem.findtext('target_exons_sam')
        print self.sam_target_exons
        if self.sam_target_exons=="":
            self.sam_target_exons=None


        # gatk params
        elem = root.find('gatk')
        self.gatk_path = elem.get('path')
        self.gatk_resource_path = elem.findtext('resource_path')
        self.gatk_use_mem = int(elem.findtext('use_mem'))
        self.gatk_num_cores = int(elem.findtext('num_cores'))
        self.gatk_snps_fdr_threshold = float(elem.findtext('snps_fdr_threshold'))
        self.gatk_recal_anal_dir = elem.findtext('recal_anal_dir')
        self.gatk_temp_dir = elem.findtext('temp_dir')
        
        # gatk unified genotyper parameters   
        stand_call_conf=elem.findtext('stand_call_conf')#default=str(50))
        stand_emit_conf=elem.findtext('stand_emit_conf')#, default=str(10))
        dcov=elem.findtext('downsample_to_coverage')#, default=str(1000))
        trig_call_conf = elem.findtext('trig_call_conf')#, default=str(45))
        trig_emit_conf = elem.findtext('trig_emit_conf') # default=str(30))
        base_quality = elem.findtext('base_quality')# default=str(30))
        read_quality = elem.findtext('read_quality')# default=str(30))
        self.gatk_unifiedGenotyper_params=[stand_call_conf,stand_emit_conf, dcov,trig_call_conf,
                                           base_quality,read_quality,trig_emit_conf]
        

        # This field need to change and be a field at the sample level
        self.gatk_target_exons = elem.findtext('target_exons_gatk')
        if self.gatk_target_exons=="":
            self.gatk_target_exons=None
        
        # This is for realigment around known indels
        self.gatk_known_realignment_sites = elem.findtext('gatk_known_realignment_sites')
        if self.gatk_known_realignment_sites=="":
            self.gatk_known_realignment_sites=None 
        
        # For Doing analysis of the Quality recalibration Procedure by analyzing the covariates.
        self.gatk_covariates_analysis=elem.findtext('gatk_covariates_analysis')
        if self.gatk_covariates_analysis=="":
            self.gatk_covariates_analysis=False
        else:
            print self.gatk_covariates_analysis 
            self.gatk_covariates_analysis=eval(str(self.gatk_covariates_analysis))
        
        # snvs Annotation using vcfCoding tool
        elem = root.find('vcfAnnotation')
        self.vcf_annot_path = elem.get('path')
        self.annot_genelist = elem.findtext('genelist')
        
        
        # VarScan Parameters
        elem = root.find('varscan')
        if elem is not None:            
            self.varscan_path = elem.get('path')
        
            vs_min_coverage=elem.findtext('min_coverage')
            vs_min_reads_alt=elem.findtext('min_reads_alt')
            min_avgbase_quality=elem.findtext('min_avgbase_quality')
            min_var_freq=elem.findtext('min_var_freq')
            min_pvalue= elem.findtext('min_pvalue')
            self.varscan_params= [vs_min_coverage,vs_min_reads_alt,min_avgbase_quality,min_var_freq,min_pvalue]
            
        # species specific configuration 
        # TODO: simple tuple for now.. if becomes more 
        # complicated then upgrade to a class
        '''
        Genome = collections.namedtuple('Genome', ['species',
                                                   'ucsc_db',
                                                   'bwa_ref_genome',
                                                   'gatk_ref_genome',
                                                   'snpdb', 
                                                   'indeldb'])
        
        '''
        
        self.genomes = {}
        for elem in root.findall("species"):
            species = elem.get("name")
            ucsc_db = elem.findtext('ucsc_db')
            bwa_ref_genome = elem.findtext('bwa_ref_genome')
            gatk_ref_genome = elem.findtext('gatk_ref_genome')
            vcf_ref_genome = elem.findtext('vcf_ref_genome')
            snpdb = elem.findtext('snpdb')
            indeldb = elem.findtext('indeldb')
            OneKgenomes = elem.findtext('OneKgenomes')
            hapmap=elem.findtext('hapmap') 
            self.genomes[species] = type('Genome', (object,),
                                           dict(species=species,
                                           ucsc_db=ucsc_db,
                                           bwa_ref_genome=bwa_ref_genome,
                                           gatk_ref_genome=gatk_ref_genome,
                                           vcf_ref_genome=vcf_ref_genome,
                                           snpdb=snpdb,
                                           indeldb=indeldb,
                                           OneKgenomes=OneKgenomes,
                                           hapmap=hapmap))


class LaneConfig(object):
    '''
    Base class to represent information about individual sequencing
    machine units, or "lanes".
    
    The structure includes job, sample, and sequencing metadata sections
    as well as a 'fastq_files' section that points to the location of the
    raw sequences to be analyzed
    
    Specific applications should inherit the base class and add
    application-specific attributes
    '''
    # Hard-coded file names (establishes naming conventions)
    FASTQC_DATA_FILE = 'fastqc_data.txt'
    FASTQC_REPORT_FILE = 'fastqc_report.html'

    def from_xml(self, xmlfile, root_dir=None):
        tree = etree.parse(xmlfile)
        root = tree.getroot()
        self.from_xml_elem(root, root_dir)
        
    def from_xml_elem(self, xmlelem, root_dir=None):
        '''
        parses an XML element tree and sets up parameter
        attributes
        
        xmlelem: root element in tree
        root_dir: (optional) path to root_dir on the local
        machine running the pipeline.  this is added as a
        prefix to all pipeline output paths
        '''        
        # job name
        self.name = xmlelem.get("name")
        # base output dir
        if root_dir is not None:
            if not os.path.exists(root_dir):
                logging.error("Cannot create %s because 'root_dir' %s does not exist" % 
                              (self.__class__.__name__, root_dir))
                return
            self.output_dir = os.path.join(root_dir, self.name)
        else:
            self.output_dir = self.name
        
        # final output destination (not local working dir)
        output_dir_elem = xmlelem.find("output_dir")
        self.dst_output_dir = output_dir_elem.text
        self.remote = parse_bool(output_dir_elem.get("remote"))
        self.remote_ip = None
        if self.remote:
            self.remote_ip = output_dir_elem.get("ip")
        # parse sample parameters
        self.run_id = xmlelem.findtext("run_id")
        self.lib_id = xmlelem.findtext("lib_id")
        self.sample_id = xmlelem.findtext("sample_id")
        self.species = xmlelem.findtext("species")
        self.sample_description = xmlelem.findtext("description")        
        # parse sequencing parameters
        self.seq_center = xmlelem.findtext("center")
        self.seq_date = xmlelem.findtext("date")
        self.seq_platform = xmlelem.findtext("platform")
        
        if self.seq_platform == 'illumina_hiseq' or self.seq_platform == 'illumina' or \
           self.seq_platform =='gaII':
            self.seq_platform = 'ILLUMINA'
        
        self.read_length = int(xmlelem.findtext("read_length"))
        self.fragment_layout = xmlelem.findtext("fragment_layout")
        self.fragment_length_mean = int(xmlelem.findtext("fragment_length_mean"))
        self.quality_scores = xmlelem.findtext("quality_scores")
        # sequence data files
        # parse fastq files
        self.output_lanedir = os.path.join(root_dir, "lane/"+self.name)
        
        src_fastq_files = {}
        dst_fastq_files = {}
        fastq_files = {}
        fastqc_data = {}
        fastqc_html = {}
        for mate_elem in xmlelem.findall("fastq_files/*"):
            mate = int(mate_elem.get("mate"))
            src_fastq_files[mate] = mate_elem.text
            # Particular to BWA: The mate file could be .gz or fq. and even .fq but really .gz when as symbolic link. This ned to be fix according to aligner
            # IN order to allow for fq o gz use without uncompressin the data
            #ext = os.path.splitext(mate_elem.text)[1]
            ext = '.fq'
            # Note 1-17-2011 alebalbin: modification of output_dir by output_lanedir
            # the fastqc files were not generated because outdir does not exist

            dst_fastq_files[mate] = os.path.join(self.output_dir, "mate%d%s" % (mate,ext))
            fastq_files[mate] = os.path.join(self.output_dir, "mate%d.fq" % (mate))
            fastqc_subdir = os.path.join(self.output_dir, '%s_fastqc' % os.path.basename(fastq_files[mate]))
            fastqc_data[mate] = os.path.join(fastqc_subdir, self.FASTQC_DATA_FILE)
            fastqc_html[mate] = os.path.join(fastqc_subdir, self.FASTQC_REPORT_FILE) 
        self.src_fastq_files = [src_fastq_files[mate] for mate in xrange(len(src_fastq_files))]
        self.dst_fastq_files = [dst_fastq_files[mate] for mate in xrange(len(dst_fastq_files))]
        self.fastq_files = [fastq_files[mate] for mate in xrange(len(fastq_files))]
        self.fastqc_data_files = [fastqc_data[mate] for mate in xrange(len(fastqc_data))]
        self.fastqc_html_files = [fastqc_html[mate] for mate in xrange(len(fastqc_data))]      

    def to_xml(self):
        root = etree.Element("lane")
        root.set("name", self.name)
        elem = etree.SubElement(root, "output_dir")
        elem.text = self.dst_output_dir
        elem.set("remote", str(self.remote))
        if self.remote:
            elem.set("ip", self.remote_ip)        
        # sample parameters
        elem = etree.SubElement(root, "run_id")
        elem.text = self.run_id
        elem = etree.SubElement(root, "lib_id")
        elem.text = self.lib_id 
        elem = etree.SubElement(root, "sample_id")
        elem.text = self.sample_id
        elem = etree.SubElement(root, "species")    
        elem.text = self.species
        elem = etree.SubElement(root, "description")
        elem.text = self.sample_description
        # sequence data parameters
        elem = etree.SubElement(root, "center")
        elem.text = self.seq_center
        elem = etree.SubElement(root, "date")
        elem.text = self.seq_date
        elem = etree.SubElement(root, "platform")
        elem.text = self.seq_platform
        elem = etree.SubElement(root, "read_length")
        elem.text = str(self.read_length)
        elem = etree.SubElement(root, "fragment_layout")
        elem.text = self.fragment_layout
        elem = etree.SubElement(root, "fragment_length_mean")
        elem.text = str(self.fragment_length_mean)
        elem = etree.SubElement(root, "quality_scores")
        elem.text = self.quality_scores
        # sequence data files
        elem = etree.SubElement(root, "fastq_files")
        for mate,fastq_file in enumerate(self.src_fastq_files):
            file_elem = etree.SubElement(elem, "file")
            file_elem.set("mate", str(mate))
            file_elem.text = fastq_file
        # write
        return etree.tostring(root)

    def from_cfg(self, cfg_file, root_dir=None):
        '''
        parses a ConfigParser file and sets up parameter
        attributes
        
        cfg_file: file to parse
        root_dir: (optional) path to root_dir on the local
        machine running the pipeline.  this is added as a
        prefix to all pipeline output paths
        '''
        config = ConfigParser.SafeConfigParser()
        config.read(cfg_file)
        # parse job parameters
        self.name = config.get("job", "name")
        # base output dir
        if root_dir is not None:
            self.output_dir = os.path.join(root_dir, self.name)
        else:
            self.output_dir = self.name
        # final output destination (not local working dir)
        self.dst_output_dir = config.get("job", "outputdir")
        self.remote = config.getboolean("job", "remote")
        self.remote_ip = None
        if config.has_option("job", "remote_ip"):            
            self.remote_ip = config.get("job", "remote_ip")
        # parse sample parameters
        self.run_id = config.get("sample", "run_id")
        self.lib_id = config.get("sample", "lib_id")
        self.sample_id = config.get("sample", "sample_id")
        self.species = config.get("sample", "species")
        self.sample_description = config.get("sample", "description")        
        # parse sequencing parameters
        self.seq_center = config.get("sequencing", "center")
        self.seq_date = config.get("sequencing", "date")
        self.seq_platform = config.get("sequencing", "platform")
        self.read_length = config.getint("sequencing", "read_length")
        self.fragment_layout = config.get("sequencing", "fragment_layout")
        self.fragment_length_mean = config.get("sequencing", "fragment_length_mean")
        self.quality_scores = config.get("sequencing", "quality_scores")

        fastq_list = config.items("fastq_files")
        src_fastq_files = {}
        dst_fastq_files = {}
        fastq_files = {}
        fastqc_data = {}
        fastqc_html = {}
        for mate_string,fastq_file in fastq_list:
            m = re.match(r'mate(\d)', mate_string)
            mate_num = int(m.group(1))
            src_fastq_files[mate_num] = fastq_file
            ext = os.path.splitext(fastq_file)[1]
            dst_fastq_files[mate_num] = os.path.join(self.output_dir, 'mate%d%s' % (mate_num,ext))
            fastq_files[mate_num] = os.path.join(self.output_dir, "mate%d.fq" % (mate_num))            
            fastqc_subdir = os.path.join(self.output_dir, '%s_fastqc' % os.path.basename(fastq_files[mate_num]))
            fastqc_data[mate_num] = os.path.join(fastqc_subdir, self.FASTQC_DATA_FILE)
            fastqc_html[mate_num] = os.path.join(fastqc_subdir, self.FASTQC_REPORT_FILE) 
        self.src_fastq_files = [src_fastq_files[mate] for mate in xrange(len(src_fastq_files))]
        self.dst_fastq_files = [dst_fastq_files[mate] for mate in xrange(len(dst_fastq_files))]
        self.fastq_files = [fastq_files[mate] for mate in xrange(len(fastq_files))]
        self.fastqc_data_files = [fastqc_data[mate] for mate in xrange(len(fastqc_data))]
        self.fastqc_html_files = [fastqc_html[mate] for mate in xrange(len(fastqc_data))]


class ExomeLaneConfig(LaneConfig):
    '''
    Supplements the LaneConfig object with Exome information
    '''
    def from_xml_elem(self, xmlelem, root_dir=None):
        # parent class init
        super(ExomeLaneConfig, self).from_xml_elem(xmlelem, root_dir)
        # bwa alignment dir
        self.align_dir = os.path.join(self.output_dir, "align")
        # bwa output SAM file
        self.align_sam_file = os.path.join(self.align_dir, "aligned_reads.sam")
        self.align_bam_file = os.path.join(self.align_dir, "aligned_reads.bam")

        # gatk output dir
        self.gatk_dir = os.path.join(self.output_dir, "gatk")
        #self.gatk_recal_dir = os.path.join(self.gatk_dir, "recal_analysis")
        
        #Specific files related to gatk preprocessing
        self.realigned_bam_file = os.path.join(self.gatk_dir,"realigned_reads.bam")
        self.matefixed_bam_file = os.path.join(self.gatk_dir,"matefixed.bam")
        self.markdup_bam_file = os.path.join(self.gatk_dir,"markdup.bam")
        ## Base quality files at the lane level
        self.recal_anal_dir = os.path.join(self.gatk_dir,'recal_anal_dir')
        self.pre_recal_table_file = os.path.join(self.gatk_dir,"pre_recal_table.csv")
        self.recal_bam_file_name = os.path.join(self.gatk_dir,"recalibrated.bam")
        self.recal_bam_indexed_name = os.path.join(self.gatk_dir,"recalibrated_sorted.bam")
        self.post_recal_table_file = os.path.join(self.gatk_dir,"post_recal_table.csv")

        

class ExomeSampleConfig(object):
    '''
    Sample configuration information
    
    LaneConfig objects are embedded here
    '''
    def from_xml_elem(self, xmlelem, sample_dir=None, lane_dir=None):
        self.name = xmlelem.get("name")
        self.patient_id = xmlelem.findtext("patient_id")
        self.category = xmlelem.findtext("category")
        # base sample output dir
        if sample_dir is not None:
            self.output_dir = os.path.join(sample_dir, self.name)
        else:
            self.output_dir = self.name
        # lanes
        self.lanes = []
        for lane_elem in xmlelem.findall("lane"):
            lane_config = ExomeLaneConfig()
            lane_config.from_xml_elem(lane_elem, lane_dir)            
            self.lanes.append(lane_config)
        
        # gatk pre-processing sample level files    
        #self.recal_anal_dir = os.path.join(self.output_dir,'recal_anal_dir')
        
        # gatk preprocessing files
        self.merged_recalibrated_bam=os.path.join(self.output_dir,self.name+"_merged.recal.bam")
        self.merged_markdup_file=os.path.join(self.output_dir,self.name+"_merged.recal.mdup.bam")
        self.sorted_mmarkdup_bam=os.path.join(self.output_dir,self.name+"_merged.recal.smdup.bam")
        # Temp 12-23-10 delete and leave the one above
        #self.sorted_mmarkdup_bam=os.path.join(self.output_dir,self.name+"_merged_recal.markdup.psorted.bam")
        # quick pre-processing in order to call with samtools
        self.merged_raw_bam=os.path.join(self.output_dir,self.name+"_merged.quick.bam")
        self.merged_mmdup_raw_bam=os.path.join(self.output_dir,self.name+"_merged.quick.mdup.bam")
        self.sorted_quickmmdup_bam=os.path.join(self.output_dir,self.name+"_merged.quick.smdup.bam")
        
        # Traditional mpileup files
        self.mpileup_file=os.path.join(self.output_dir,self.name+".mpileup")
        self.benign_reg_mpileup_file=os.path.join(self.output_dir,self.name+'_benign.pileup')
        self.tumor_reg_mpileup_file=os.path.join(self.output_dir,self.name+'_tumor.pileup')
        self.benign_bcf_mpileup_file=os.path.join(self.output_dir,self.name+'_benign.bcf')
        self.tumor_bcf_mpileup_file=os.path.join(self.output_dir,self.name+'_tumor.bcf')
        

class ExomeAnalysisConfig(object):
    '''
    Data structure to store analysis, sample, and lane
    configurations for Exome analysis pipelines
    '''
    def not_quick(self,use_recal):
        if use_recal:
            sam_calls_dir = os.path.join(self.analysis_dir,"sam_calls")
        else:
            sam_calls_dir = os.path.join(self.analysis_dir,"sam_calls_quick")
        
        return sam_calls_dir
        
    def from_xml(self, xmlfile, root_dir=None):
        # parse the XML
        tree = etree.parse(xmlfile)
        root = tree.getroot()
        self.name = root.get("name") 

        if root_dir is None:
            root_dir = root.find("remote_output_dir")
            if root_dir is not None:
                self.remote_root_dir = root_dir.get('path')
                root_dir = self.remote_root_dir
                print root_dir
            else:
                self.remote_root_dir = os.getcwd()
                root_dir = os.getcwd()
        
        remote_root_dir = root.find("remote_output_dir")
        if remote_root_dir is not None:
            self.remote_root_dir = remote_root_dir.get('path')
            
        # setup analysis, sample, and lane root directories
        analysis_dir = os.path.join(root_dir, "analysis")
        sample_dir = os.path.join(root_dir, "sample")
        lane_dir = os.path.join(root_dir, "lane")

        # base analysis output dir
        # All this folders should be created before the analysis starts
        self.analysis_dir=analysis_dir
        self.sample_dir=sample_dir
        self.lane_dir=lane_dir
        self.output_dir = os.path.join(analysis_dir, self.name)
        # To realign only at known indel sites
        self.known_realignment_sites =  os.path.join(self.output_dir, self.name+"_knownSitesForRealigning.intervals")  
        
        ####  mpileup files: Traditional ones for coverage and with genotypes for sam_tools
        self.mpileup_file=os.path.join(self.output_dir,self.name+".mpileup")
        self.benign_reg_mpileup_file=os.path.join(self.output_dir,self.name+'_benign.pileup')
        self.tumor_reg_mpileup_file=os.path.join(self.output_dir,self.name+'_tumor.pileup')

        ###########
        # final output destination (not local working dir)
        output_dir_elem = root.find("output_dir")
        self.dst_output_dir = output_dir_elem.text
        self.remote = parse_bool(output_dir_elem.get("remote"))
        self.remote_ip = None
        if self.remote:
            self.remote_ip = output_dir_elem.get("ip")
        # samples
        self.samples = []
        for sample_elem in root.findall("sample"):
            sample_config = ExomeSampleConfig()
            sample_config.from_xml_elem(sample_elem, sample_dir, lane_dir)
            self.samples.append(sample_config)
        # check species
        species = set()
        for sample in self.samples:
            for lane in sample.lanes:
                species.add(lane.species)
        if len(species) > 1:
            logging.error("Analysis file has data from multiple species")
            assert False
        self.species = list(species)[0]
        
        ##################GATK calls files #####################
        self.gatk_calls_dir = os.path.join(self.analysis_dir,"gatk_calls")
        self.snps_raw_vcf=os.path.join(self.gatk_calls_dir,self.name+'.gatk_snps.raw.vcf')
  
        # Calling raw indels and creating a mask indel file
        self.indels_vcf = os.path.join(self.gatk_calls_dir,self.name+'.gatk_indels.raw.vcf')
        self.indels_bed =os.path.join(self.gatk_calls_dir,self.name+'.gatk_indels.raw.bed')
        self.indels_masked_bed =os.path.join(self.gatk_calls_dir,self.name+'.gatk_indels.mask.bed')
        
        #### Multisample analysis ######
        self.snps_basic_vcf =os.path.join(self.gatk_calls_dir,self.name+'.gatk_snps.filtered.basic.vcf')
        # Method I)  Hard Filtering
        # Generates snps_hard_vcf file
        self.snps_hard_vcf =os.path.join(self.gatk_calls_dir,self.name+'.gatk_snps.filtered.hard.vcf')
        self.snps_hard_annotated_vcf = os.path.join(self.gatk_calls_dir,self.name+
                                                    '.gatk_snps.filtered.hard_annotated.vcf')
        self.snps_hard_annotated_log = os.path.join(self.gatk_calls_dir,self.name+
                                                    '.gatk_snps.filtered.hard_annotated.log')
        
        # Method II) Model based calling
        # Three steps: a) cluster generation, b) calculation of variant posterior probability
        # c) select variants according to specific fdr
        # Generates snps_model_vcf file
        self.cluster_ouput_file =os.path.join(self.gatk_calls_dir,self.name+'.gatk_snps.model.clusters')
        self.snps_variantRecal_vcf =os.path.join(self.gatk_calls_dir,self.name+'.gatk_snps.filtered.variantRecal.vcf')
        self.snps_model_vcf = os.path.join(self.gatk_calls_dir,self.name+'.gatk_snps.filtered.model.vcf')
        self.tranches_file =os.path.join(self.gatk_calls_dir,self.name+'.gatk_snps.filtered.model.tranches')
        
        self.snps_model_annotated_vcf=os.path.join(self.gatk_calls_dir,self.name+
                                                   '.gatk_snps.filtered.model_annotated.vcf')
        self.snps_model_annotated_log=os.path.join(self.gatk_calls_dir,self.name+
                                                   '.gatk_snps.filtered.model_annotated.log')

        
        #### Paired-Sample analysis ######
        ### Hard Filtering
        ## Benign
        self.benign_snps_raw_vcf=os.path.join(self.gatk_calls_dir,self.name+'.gatk_snps.b.raw.vcf')
        self.benign_snps_basic_vcf =os.path.join(self.gatk_calls_dir,self.name+'.gatk_snps.b.filtered.basic.vcf')
        self.benign_snps_hard_vcf =os.path.join(self.gatk_calls_dir,self.name+'.gatk_snps.b.filtered.hard.vcf')
        self.benign_snps_hard_annotated_vcf = os.path.join(self.gatk_calls_dir,self.name+
                                                    '.gatk_snps.b.filtered.hard_annotated.vcf')
        self.benign_snps_hard_annotated_log = os.path.join(self.gatk_calls_dir,self.name+
                                                    '.gatk_snps.b.filtered.hard_annotated.log')
        
        
        ## Tumors
        self.tumor_snps_raw_vcf=os.path.join(self.gatk_calls_dir,self.name+'.gatk_snps.t.raw.vcf')
        self.tumor_snps_basic_vcf =os.path.join(self.gatk_calls_dir,self.name+'.gatk_snps.t.filtered.basic.vcf')
        self.tumor_snps_hard_vcf =os.path.join(self.gatk_calls_dir,self.name+'.gatk_snps.t.filtered.hard.vcf')
        
        self.tumor_snps_hard_annotated_vcf=os.path.join(self.gatk_calls_dir,self.name+
                                                   '.gatk_snps.t.filtered.hard_annotated.vcf')
        self.tumor_snps_hard_annotated_log=os.path.join(self.gatk_calls_dir,self.name+
                                                   '.gatk_snps.t.filtered.hard_annotated.log')
  
        # Variant Evaluation
        self.hard_merged_tumor_benign_vcf = os.path.join(self.gatk_calls_dir,self.name+'gatk_snps.m.fhard.vcf')
        self.hard_variant_eval_file_sufix = os.path.join(self.gatk_calls_dir,self.name+'gatk_snps.m.fhard.eval')
        self.hard_interesting_snps_vcf = os.path.join(self.gatk_calls_dir,self.name+'gatk_snps.m.fhard.interest.vcf')

        ### Model based
        ## Benign
        self.benign_cluster_ouput_file =os.path.join(self.gatk_calls_dir,self.name+'.gatk_snps.b.model.clusters')
        self.benign_snps_variantRecal_vcf =os.path.join(self.gatk_calls_dir,self.name+'.gatk_snps.b.filtered.variantRecal.vcf')
        self.benign_snps_model_vcf = os.path.join(self.gatk_calls_dir,self.name+'.gatk_snps.b.filtered.model.vcf')
        self.benign_tranches_file =os.path.join(self.gatk_calls_dir,self.name+'.gatk_snps.b.filtered.model.tranches')
        ### Tumors
        self.tumor_cluster_ouput_file =os.path.join(self.gatk_calls_dir,self.name+'.gatk_snps.t.model.clusters')
        self.tumor_snps_variantRecal_vcf =os.path.join(self.gatk_calls_dir,self.name+'.gatk_snps.t.filtered.variantRecal.vcf')
        self.tumor_snps_model_vcf = os.path.join(self.gatk_calls_dir,self.name+'.gatk_snps.t.filtered.model.vcf')
        self.tumor_tranches_file =os.path.join(self.gatk_calls_dir,self.name+'.gatk_snps.t.filtered.model.tranches')
        
        
        
        # Variant Evaluation
        self.model_merged_tumor_benign_vcf = os.path.join(self.gatk_calls_dir,self.name+'gatk_snps.m.fmodel.vcf')
        self.model_variant_eval_file_sufix = os.path.join(self.gatk_calls_dir,self.name+'gatk_snps.m.fmodel.eval')
        self.model_interesting_snps_vcf = os.path.join(self.gatk_calls_dir,self.name+'gatk_snps.m.fmodel.interest.vcf')
        
        self.benign_snps_model_annotated_vcf=os.path.join(self.gatk_calls_dir,self.name+
                                                   '.gatk_snps.b.filtered.model_annotated.vcf')
        self.benign_snps_model_annotated_log=os.path.join(self.gatk_calls_dir,self.name+
                                                   '.gatk_snps.b.filtered.model_annotated.log')
        
        self.tumor_snps_model_annotated_vcf=os.path.join(self.gatk_calls_dir,self.name+
                                                   '.gatk_snps.t.filtered.model_annotated.vcf')
        self.tumor_snps_model_annotated_log=os.path.join(self.gatk_calls_dir,self.name+
                                                   '.gatk_snps.t.filtered.model_annotated.log')
        
        
    def sam_tools_files(self,use_recal):
        self.sam_calls_dir = self.not_quick(use_recal)
        ###################SAMTOOLS calls files############################
        # Benign
        #self.sam_calls_dir = is_analysis_quick(self,quick)
        #self.sam_calls_dir = os.path.join(self.analysis_dir,"sam_calls")
        self.benign_bcf_mpileup_file = os.path.join(self.sam_calls_dir,self.name+'benign.sam_snps.mpileup.bcf')
        self.benign_mpileup_file = os.path.join(self.sam_calls_dir,self.name+'benign.sam_snps.mpileup.vcf')
        self.benign_filtered_snps = os.path.join(self.sam_calls_dir,self.name+'benign.sam_snps.filtered.recode.vcf')
        self.benign_snps_stats = os.path.join(self.sam_calls_dir,self.name+'benign.sam_snps.filtered.stats')
        #Tumor
        self.tumor_bcf_mpileup_file=os.path.join(self.sam_calls_dir,self.name+'tumor.sam_snps.mpileup.bcf')
        self.tumor_mpileup_file = os.path.join(self.sam_calls_dir,self.name+'tumor.sam_snps.mpileup.vcf')
        self.tumor_filtered_snps = os.path.join(self.sam_calls_dir,self.name+'tumor.sam_snps.filtered.recode.vcf')
        self.tumor_snps_stats = os.path.join(self.sam_calls_dir,self.name+'tumor.sam_snps.filtered.stats')
    
    def varscan_files(self):
        '''
        files for the varscan analysis
        '''
        self.varscan_calls_dir = os.path.join(self.analysis_dir,"varscan_calls")
        ###########
        ### VarScan
        self.snps_varscan_somatic = os.path.join(self.varscan_calls_dir,self.name+'_snps_varscan_somatic')
        self.indels_varscan_somatic = os.path.join(self.varscan_calls_dir,self.name+'_indels_varscan_somatic')        
        self.snps_varscan_single = os.path.join(self.varscan_calls_dir,self.name+'_snps_varscan_single')
        self.indels_varscan_single = os.path.join(self.varscan_calls_dir,self.name+'_indels_varscan_single')




if __name__ == '__main__':
    import sys
    root_dir = sys.argv[2]
    analysis = ExomeAnalysisConfig()
    analysis.from_xml(sys.argv[1], root_dir)
    print analysis.name
    print analysis.output_dir
    for sample in analysis.samples:
        print sample.name
        print sample.patient_id
        print sample.category
        print sample.output_dir
        for lane in sample.lanes:
            print lane.name
            print lane.output_dir
            print lane.align_bam_file

